ESS 575 Models for Ecological Data
Swiss Birds Occupancy Modeling
April 26, 2017
I. Motivation
Modeling presence or absence is a classic problem involving mixture models, specifically random variables that are zero-inflated. Extra zeros are encountered when we model presence or absence because zeros arise from two conditions: truly absent individuals and individuals present but undetected. This means we need a model for the process that controls occupancy, the true state, and model of the data that accounts for detection. This is often our starting point in Bayesian analysis – there is a true, unobserved state we seek to understand using a model of a process. We take imperfect observations on that state and must correct them using a model of the data.
II. Problem
A fundamental question in landscape ecology seeks to understand how landscape structure shapes variation in the abundance of species. We will use data from the Swiss Survey of Common Breeding Birds, courtesy of Royle and Dorazio (2008), to model habitat occupancy by a common, resident bird in the Swiss Alps, the willow tit (Parus montanus). The data come from annual surveys of one km2 quadrats distributed across Switzerland (Fig. 1). Surveys are conducted during the breeding season on three separate days, but some quadrats have missing data so that the number of replicate observations is fewer than three.
During each survey, an observer records every visual or acoustic detection of a breeding species (we do not differentiate between these two types of detection in this problem) and marks its location using a global positioning system or, in earlier years, a paper map. We assume that the true state (occupied or unoccupied) does not change among sample dates, an assumption known as closure. This assumption is reasonable because we are observing a resident species during the breeding season.
Fig. 1. The willow tit (left, c/o of Francis C. Franklin) is one of 70 bird species that are surveyed annually for abundance in 267 1 km2 sampling units distributed across Switzerland (right, c/o the Swiss Ornithological Institute).
We want to understand the influence of forest cover and elevation on the distribution of the willow tit. The data frame SwissBirds has the number of times a quadrat (quadrat) was searched (numberVisits) and the number of times willow tits were detected (numberDetections). We have covariates on forest canopy cover (forestCover) as well as elevation in meters (elevation) for each quadrat surveyed. Develop a model of the influence of forest cover and elevation on the distribution of willow tits. Your model should allow estimation of the optimum elevation of willow tit habitat at the mean forest cover, where optimum elevation is defined as the elevation where probability of occupancy is maximum.
- Diagram the network of knowns and unknowns.
- Write out a mathematical expression for the posterior and the joint distribution.
- Approximate the marginal posterior distributions of the model parameters using JAGS. Some hints:
- You will need to standardize the covariates by subtracting the mean and dividing by the standard deviation for each observation in the elevation and forest cover data. Use the
scalefunction to scale these covariates (this will drastically speed convergence). - You must give initial values of 1 to all unknown 0 or 1 z states.
{ # Extra bracket needed only for R markdown files
sink("SwissBirds.R")
cat("
model{
# priors
p ~ dbeta(1, 1)
for (i in 1:4){
beta[i] ~ dnorm(0, .368)
}
# likelihood
for (i in 1:N){
z[i] ~ dbern(psi[i])
logit(psi[i]) <- beta[1] + beta[2]*elevation[i] + beta[3]*elevation[i]^2 + beta[4]*forestCover[i]
y[i] ~ dbin(z[i] * p, n[i])
}
# derived quantities
elevationMaxScaled <- -beta[2]/(2 * beta[3])
elevationMax <- elevationMaxScaled * sdElevation + muElevation
for (j in 1:length(elevationPredict)){
#note that the x's are standardized on the R side
logit(psiPredict[j]) <- beta[1] + beta[2]*elevationPredict[j] + beta[3]*elevationPredict[j]^2
}
}
", fill = TRUE)
sink()
} # Extra bracket needed only for R markdown fileslibrary(ESS575)
library(rjags)## Loading required package: coda## Linked to JAGS 4.2.0## Loaded modules: basemod,bugsset.seed(1)
elevation <- scale(SwissBirds$elevation)
muElevation <- mean(SwissBirds$elevation)
sdElevation <- sd(SwissBirds$elevation)
elevation_Predict = seq(500,2500,10) # for plotting probability of occupancy vs unstandardized elevation
elevationPredict_stand = (elevation_Predict - muElevation)/sdElevation #Standardized for making predictions of probability of occupancy at new elevations
forestCover <- scale(SwissBirds$forestCover)
muforestCover <- mean(SwissBirds$forestCover)
sdforestCover <- sd(SwissBirds$forestCover)
inits = list(
list(z = rep(1, nrow(SwissBirds)), p = runif(1, 0, 1), beta = runif(4, -2, 2)),
list(z = rep(1, nrow(SwissBirds)), p = runif(1, 0, 1), beta = runif(4, -2, 2)),
list(z = rep(1, nrow(SwissBirds)), p = runif(1, 0, 1), beta = runif(4, -2, 2)))
data = list(
N = nrow(SwissBirds),
elevation = as.double(elevation),
forestCover = as.double(forestCover),
n = as.double(SwissBirds$numberVisits),
y = as.double(SwissBirds$numberDetections),
muElevation = as.double(muElevation),
sdElevation = as.double(sdElevation),
elevationPredict = as.double(elevationPredict_stand)) #note standardized for predictions
n.adapt = 5000
n.update = 10000
n.iter = 10000
jm = jags.model("SwissBirds.R", data = data, inits = inits, n.chains = length(inits), n.adapt = n.adapt)## Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph information:
## Observed stochastic nodes: 237
## Unobserved stochastic nodes: 242
## Total graph size: 4128
##
## Initializing modelupdate(jm, n.iter = n.update)
zm = coda.samples(jm, variable.names = c("p", "beta"), n.iter = n.iter, n.thin = 1)
zj = jags.samples(jm, variable.names = c("elevationMax", "psiPredict"), n.iter = n.iter, n.thin = 1)- Summarize the parameters and check chains for convergence.
summary(zm)##
## Iterations = 15001:25000
## Thinning interval = 1
## Number of chains = 3
## Sample size per chain = 10000
##
## 1. Empirical mean and standard deviation for each variable,
## plus standard error of the mean:
##
## Mean SD Naive SE Time-series SE
## beta[1] -0.2003 0.27509 0.0015882 0.0042699
## beta[2] 1.9680 0.29755 0.0017179 0.0042504
## beta[3] -1.0961 0.26586 0.0015349 0.0048402
## beta[4] 0.8435 0.23296 0.0013450 0.0029456
## p 0.7873 0.02873 0.0001659 0.0002347
##
## 2. Quantiles for each variable:
##
## 2.5% 25% 50% 75% 97.5%
## beta[1] -0.7425 -0.3844 -0.2000 -0.01667 0.3400
## beta[2] 1.4226 1.7598 1.9562 2.16200 2.5764
## beta[3] -1.6428 -1.2699 -1.0852 -0.91232 -0.6026
## beta[4] 0.4070 0.6817 0.8389 0.99857 1.3152
## p 0.7283 0.7684 0.7886 0.80731 0.8402plot(zm)
gelman.diag(zm)## Potential scale reduction factors:
##
## Point est. Upper C.I.
## beta[1] 1 1
## beta[2] 1 1
## beta[3] 1 1
## beta[4] 1 1
## p 1 1
##
## Multivariate psrf
##
## 1heidel.diag(zm)## [[1]]
##
## Stationarity start p-value
## test iteration
## beta[1] passed 1 0.494
## beta[2] passed 1 0.175
## beta[3] passed 1 0.186
## beta[4] passed 1 0.768
## p passed 1 0.173
##
## Halfwidth Mean Halfwidth
## test
## beta[1] passed -0.207 0.014397
## beta[2] passed 1.964 0.013477
## beta[3] passed -1.089 0.015544
## beta[4] passed 0.845 0.009996
## p passed 0.787 0.000775
##
## [[2]]
##
## Stationarity start p-value
## test iteration
## beta[1] passed 1 0.174
## beta[2] passed 1 0.686
## beta[3] passed 1 0.778
## beta[4] passed 3001 0.908
## p passed 1 0.435
##
## Halfwidth Mean Halfwidth
## test
## beta[1] passed -0.206 0.014248
## beta[2] passed 1.972 0.015549
## beta[3] passed -1.095 0.017689
## beta[4] passed 0.839 0.011342
## p passed 0.787 0.000811
##
## [[3]]
##
## Stationarity start p-value
## test iteration
## beta[1] passed 1 0.805
## beta[2] passed 1 0.557
## beta[3] passed 1 0.949
## beta[4] passed 1 0.973
## p passed 1 0.556
##
## Halfwidth Mean Halfwidth
## test
## beta[1] passed -0.188 0.014836
## beta[2] passed 1.968 0.014185
## beta[3] passed -1.105 0.015984
## beta[4] passed 0.838 0.010441
## p passed 0.787 0.000804- What can you conclude about the relative importance of elevation and forest cover in controlling the bird’s distribution? Plot the probability of occupancy as function of elevation at the mean of forest cover. Find the optimum elevation of willow tit habitat at the mean forest cover, where optimum elevation is defined as the elevation where probability of occupancy is maximum. Plot a normalized histogram of MCMC output for the optimum elevation at the average forest cover. Overlay 95% credible intervals on the optimum elevation.
We are entitled to conclude that elevation is roughly twice as important as forest cover because its coefficient the elevation coefficient is twice as large and because the data have been standardized.
You need to formulate a deterministic model that has a maximum and that is relatively easy to differentiate, something like \(\mu_i = \beta_0+ \beta_1x_{1,i}+\beta_2x_{2,i}^2+\beta_3x_{3,1}\) where \(\mathbf{x}_2\) are the data on elevation and \(\mathbf{x}_3\) are the data on forest cover, when = 0 at the mean. Dropping the i subscript for more general notation, we now have \[\mu=\beta_0+ \beta_1x_1+\beta_2x_2^2.\]
To find the maximum, we take the first derivative \(\frac {d\mu}{dx}=\beta_1+ 2\beta_2 x_{1}\), set it to zero \(0=\beta_2+ 2\beta x_1\), and solve for \(x_1\), \(x_1=\frac{-\beta_1}{2\beta_2}\). You include this result as a derived quantity in your JAGS code.
par(mfrow = c(1, 2))
psi <- summary(zj$psiPredict, quantile, c(.025, .5, .975))$stat
plot(elevation_Predict, psi[2,], type = "l", ylim = c(0,1), xlab = "Elevation (m)", lwd = 2,
ylab ="Probability of occupancy")
lines(elevation_Predict, psi[1,], type = "l", lty = "dashed", lwd = 2)
lines(elevation_Predict, psi[3,], type = "l", lty = "dashed", lwd = 2)
abline(v = mean(zj$elevationMax), lwd = 2)
text(1300, .9, "Optimum elevation", cex = .8)
hist(zj$elevationMax, breaks = 500, main = "", xlab = "Elevation (m)", xlim = c(1500, 2500),
freq = FALSE, col = "azure1")
abline(v = summary(zj$elevationMax, quantile, c(.975))$stat, lty = "dashed", lwd = 2)
abline(v = summary(zj$elevationMax, quantile, c(.025))$stat, lty = "dashed", lwd = 2)
Fig.2. Probability of occupancy at the mean forest cover as a function of elevation (left panel) and the posterior density of the optimum elevation at the mean forest cover (right panel). Dashed give are .025 and .975 quantiles, also known as .95 equal-tailed Bayesian credible intervals.
III. References
Royle, J.A., and R.M. Dorazio. 2008. Hierarchical Modeling and Inference in Ecology: The Analysis of Data from Populations, Metapopulations, and Communities. Academic Press, London, United Kingdom.